Повторение примеров

Решение обыкновенных дифференциальных уравне

Модель экспоненциального роста

using DifferentialEquations
# задаём описание модели с начальными условиями:
a = 0.98
f(u,p,t) = a*u
u0 = 1.0
# задаём интервал времени:
tspan = (0.0,1.0)
# решение:
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)
# подключаем необходимые пакеты:
using Plots
# строим графики:
plot(sol, linewidth=5,title="Модель экспоненциального роста", xaxis="Время",yaxis="u(t)",label="u(t)")
plot!(sol.t, t->1.0*exp(a*t),lw=3,ls=:dash,label="Аналитическое решение")

# задаём точность решения:
sol = solve(prob,abstol=1e-8,reltol=1e-8)
println(sol)
# строим график:
plot(sol, lw=2, color="black", title="Модель экспоненциального роста", xaxis="Время",yaxis="u(t)",label="Численное решение")
plot!(sol.t, t->1.0*exp(a*t),lw=3,ls=:dash,color="red",label="Аналитическое решение")

ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CompositeAlgorithm{Tuple{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5P{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5P{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Vern7ConstantCache, OrdinaryDiffEq.Rosenbrock5ConstantCache{SciMLBase.TimeDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, SciMLBase.UDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}, Float64, OrdinaryDiffEq.StaticWOperator{true, Float64}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5P{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}, SciMLBase.DEStats, Vector{Int64}}([1.0, 1.0412786454705882, 1.154721667167105, 1.3239023451498872, 1.536369419677937, 1.8214715431940545, 2.187114005151482, 2.662728897933159, 2.664456241933516], nothing, nothing, [0.0, 0.04127492324135852, 0.14679523890358218, 0.28630989526009537, 0.43818583001107836, 0.611882361849277, 0.7985539478479285, 0.9993382636667485, 1.0], [[1.0], [0.98, 0.9801982021814051, 0.9843259088450718, 0.9864960397068047, 0.9982235849699092, 1.0044615253487053, 1.0156761793421412, 1.0173619945010424, 1.0204530620487353, 1.0204531254788294], [1.0204530725611765, 1.0209806973598188, 1.0320082607587013, 1.03783519380158, 1.0696697047343133, 1.0868419096646311, 1.1181335167357913, 1.1228839561927564, 1.1316267628379957, 1.1316296044455856], [1.131627233823763, 1.1324008388885507, 1.1486000430640122, 1.1571826056824441, 1.2043433970666335, 1.2299724884788301, 1.2770083451466132, 1.2841859723003326, 1.2974226977952101, 1.2974323547575104], [1.2974242982468895, 1.2983898311344073, 1.3186218509251284, 1.3293515279843413, 1.3884344417206802, 1.42062931021883, 1.4798681016827064, 1.4889248128328403, 1.5056394522771077, 1.505655014644132], [1.5056420312843781, 1.5069235027995953, 1.5338082717543726, 1.5480907641278678, 1.6270290248411337, 1.6702488207053072, 1.7501372010794107, 1.762391153483508, 1.7850369857394834, 1.7850679249179284], [1.7850421123301734, 1.7866748738768703, 1.820954186120135, 1.839183776362968, 1.940159444284083, 1.9956020144031121, 2.098361267550903, 2.1141540608294553, 2.1433636124760094, 2.1434125768523966], [2.1433717250484525, 2.145480466632713, 2.1897874094236114, 2.213376094416359, 2.3443496552456504, 2.416485322840126, 2.5505783188566515, 2.571230230985217, 2.6094612750321717, 2.609540018635808], [2.609474319974496, 2.60948278121594, 2.609658593516366, 2.609750735167298, 2.6102452529283604, 2.610505965171262, 2.6109706963387507, 2.6110401192280173, 2.611167117094843, 2.6111671170948547]], ODEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}(ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}(f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, nothing), 1.0, (0.0, 1.0), SciMLBase.NullParameters(), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}(), SciMLBase.StandardODEProblem()), CompositeAlgorithm(; algs = (Vern7(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false), lazy = true,), Rodas5P(; linsolve = nothing, precs = DEFAULT_PRECS,)), choice_function = OrdinaryDiffEq.AutoSwitchCache{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5P{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}(-8, 8, Vern7(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false), lazy = true,), Rodas5P(; linsolve = nothing, precs = DEFAULT_PRECS,), false, 10, 3, 9//10, 9//10, 2, false, 5),), OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Vern7ConstantCache, OrdinaryDiffEq.Rosenbrock5ConstantCache{SciMLBase.TimeDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, SciMLBase.UDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}, Float64, OrdinaryDiffEq.StaticWOperator{true, Float64}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5P{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}(ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}(f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, nothing), [1.0, 1.0412786454705882, 1.154721667167105, 1.3239023451498872, 1.536369419677937, 1.8214715431940545, 2.187114005151482, 2.662728897933159, 2.664456241933516], [0.0, 0.04127492324135852, 0.14679523890358218, 0.28630989526009537, 0.43818583001107836, 0.611882361849277, 0.7985539478479285, 0.9993382636667485, 1.0], [[1.0], [0.98, 0.9801982021814051, 0.9843259088450718, 0.9864960397068047, 0.9982235849699092, 1.0044615253487053, 1.0156761793421412, 1.0173619945010424, 1.0204530620487353, 1.0204531254788294], [1.0204530725611765, 1.0209806973598188, 1.0320082607587013, 1.03783519380158, 1.0696697047343133, 1.0868419096646311, 1.1181335167357913, 1.1228839561927564, 1.1316267628379957, 1.1316296044455856], [1.131627233823763, 1.1324008388885507, 1.1486000430640122, 1.1571826056824441, 1.2043433970666335, 1.2299724884788301, 1.2770083451466132, 1.2841859723003326, 1.2974226977952101, 1.2974323547575104], [1.2974242982468895, 1.2983898311344073, 1.3186218509251284, 1.3293515279843413, 1.3884344417206802, 1.42062931021883, 1.4798681016827064, 1.4889248128328403, 1.5056394522771077, 1.505655014644132], [1.5056420312843781, 1.5069235027995953, 1.5338082717543726, 1.5480907641278678, 1.6270290248411337, 1.6702488207053072, 1.7501372010794107, 1.762391153483508, 1.7850369857394834, 1.7850679249179284], [1.7850421123301734, 1.7866748738768703, 1.820954186120135, 1.839183776362968, 1.940159444284083, 1.9956020144031121, 2.098361267550903, 2.1141540608294553, 2.1433636124760094, 2.1434125768523966], [2.1433717250484525, 2.145480466632713, 2.1897874094236114, 2.213376094416359, 2.3443496552456504, 2.416485322840126, 2.5505783188566515, 2.571230230985217, 2.6094612750321717, 2.609540018635808], [2.609474319974496, 2.60948278121594, 2.609658593516366, 2.609750735167298, 2.6102452529283604, 2.610505965171262, 2.6109706963387507, 2.6110401192280173, 2.611167117094843, 2.6111671170948547]], [1, 1, 1, 1, 1, 1, 1, 1, 1], true, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Vern7ConstantCache, OrdinaryDiffEq.Rosenbrock5ConstantCache{SciMLBase.TimeDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, SciMLBase.UDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}, Float64, OrdinaryDiffEq.StaticWOperator{true, Float64}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5P{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}((OrdinaryDiffEq.Vern7ConstantCache(), OrdinaryDiffEq.Rosenbrock5ConstantCache{SciMLBase.TimeDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, SciMLBase.UDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}, Float64, OrdinaryDiffEq.StaticWOperator{true, Float64}, Nothing}(SciMLBase.TimeDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}(ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}(f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, nothing), 1.0, SciMLBase.NullParameters()), SciMLBase.UDerivativeWrapper{false, ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}(ODEFunction{false, SciMLBase.AutoSpecialize, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}(f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, nothing), 0.0, SciMLBase.NullParameters()), OrdinaryDiffEq.Rodas5Tableau{Float64, Float64}(3.0, 2.849394379747939, 0.45842242204463923, -6.954028509809101, 2.489845061869568, -10.358996098473584, 2.8029986275628964, 0.5072464736228206, -0.3988312541770524, -0.04721187230404641, -7.502846399306121, 2.561846144803919, -11.627539656261098, -0.18268767659942256, 0.030198172008377946, -14.155112264123755, -17.97296035885952, -2.859693295451294, 147.12150275711716, -1.41221402718213, 71.68940251302358, 165.43517024871676, -0.4592823456491126, 42.90938336958603, -5.961986721573306, 24.854864614690072, -3.0009227002832186, 47.4931110020768, 5.5814197821558125, -0.6610691825249471, 30.91273214028599, -3.1208243349937974, 77.79954646070892, 34.28646028294783, -19.097331116725623, -28.087943162872662, 37.80277123390563, -3.2571969029072276, 112.26918849496327, 66.9347231244047, -40.06618937091002, -54.66780262877968, -9.48861652309627, 0.21193756319429014, 0.21193756319429014, -0.42387512638858027, -0.3384627126235924, 1.8046452872882734, 2.325825639765069, 0.6358126895828704, 0.4095798393397535, 0.9769306725060716, 0.4288403609558664, 25.948786856663858, -2.5579724845846235, 10.433815404888879, -2.3679251022685204, 0.524948541321073, 1.1241088310450404, 0.4272876194431874, -0.17202221070155493, -9.91568850695171, -0.9689944594115154, 3.0438037242978453, -24.495224566215796, 20.176138334709044, 15.98066361424651, -6.789040303419874, -6.710236069923372, 11.419903575922262, 2.8879645146136994, 72.92137995996029, 80.12511834622643, -52.072871366152654, -59.78993625266729, -0.15582684282751913, 4.883087185713722), 0.0, OrdinaryDiffEq.StaticWOperator{true, Float64}(0.0), nothing)), OrdinaryDiffEq.AutoSwitchCache{Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rodas5P{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}(-8, 8, Vern7(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false), lazy = true,), Rodas5P(; linsolve = nothing, precs = DEFAULT_PRECS,), false, 10, 3, 9//10, 9//10, 2, false, 5), 1)), true, 0, SciMLBase.DEStats(82, 0, 0, 0, 0, 0, 0, 0, 8, 0, 1.0), [1, 1, 1, 1, 1, 1, 1, 1, 1], SciMLBase.ReturnCode.Success)

Система Лоренца

# задаём описание модели:
function lorenz!(du,u,p,t)
σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
du[3] = u[1]*u[2] - β*u[3]
end
# задаём начальное условие:
u0 = [1.0,0.0,0.0]
# задаём знанчения параметров:
p = (10,28,8/3)
# задаём интервал времени:
tspan = (0.0,100.0)
# решение:
prob = ODEProblem(lorenz!,u0,tspan,p)
sol = solve(prob)
# строим график:
plot(sol, vars=(1,2,3), lw=2, title="Аттрактор Лоренца", xaxis="x",yaxis="y", zaxis="z",legend=false)

┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead.
│   caller = ip:0x0
└ @ Core :-1

Без интерполяции

# отключаем интерполяцию:
plot(sol,vars=(1,2,3),denseplot=false, lw=1, title="Аттрактор Лоренца", xaxis="x",yaxis="y", zaxis="z",legend=false)

Модель Лотки–Вольтерры

using ParameterizedFunctions, DifferentialEquations, Plots;
# задаём описание модели:
lv! = @ode_def LotkaVolterra begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a b c d
# задаём начальное условие:
u0 = [1.0,1.0]
# задаём знанчения параметров:
p = (1.5,1.0,3.0,1.0)
# задаём интервал времени:
tspan = (0.0,10.0)
# решение:
prob = ODEProblem(lv!,u0,tspan,p)
sol = solve(prob)
plot(sol, label = ["Жертвы" "Хищники"], color="black", ls=[:solid :dash], title="Модель Лотки - Вольтерры", xaxis="Время",yaxis="Размер популяции")

# фазовый портрет:
plot(sol,vars=(1,2), color="black", xaxis="Жертвы",yaxis="Хищники", legend=false)

Самостоятельная работа

Модель Мальтуса

Модель Мальтуса --- модель роста численности изолированной популяции, где изменение роста популяции контролируется численностью уже существующей популяции, домноженной на коэффициент $a$, который является разницей между рождаемостью и смертностью ($b-c$). Коэффициенты $b$ и $c$ было предложено выбрать самостоятельно, и я выставлю для системы значения $b=1.09$ и $c=1.134$ (что является соответственно коэффициентами рождаемости и смертности за январь-август в 2022 году в Центральном федеральном округе РФ). Изначальная численность населения ($39 433 556$ человек) также взята из статистики Росстата за 2022 год (с учётом переписи населения).

Модель Мальтуса подразумевает, что коэффициенты рождаемости и смертности не изменяются, так что если $b$ превышает $c$, численность популяции будет расти (и наоборот).

function Maltus!(du,u,p,t)
    du[1] = (p[1]-p[2])*u[1]
end
u0 = [39433556.0]
tspan = (0.0,100.0)
p = Float64[1.09, 1.134]
prob = ODEProblem(Maltus!,u0,tspan,p)
sol = solve(prob,abstol=1e-6,reltol=1e-6, saveat=1.0)
R1 = [tu[1] for tu in sol.u]
plot(sol.t, R1, title="Модель Мальтуса", xaxis="Время",yaxis="Численность",label="u(t)", c=:blue)

anim = @animate for i in 1:length(sol.t)
    plot(sol.t[1:i], R1[1:i], title="Модель Мальтуса", xaxis="Время",yaxis="Численность",label="u(t)", c=:blue)
end
gif(anim, "presentation//image//1.gif")

[ Info: Saved animation to D:\Education\КомпПрактикумПоСтатМоделированию\labs\gitrepo\lab6\presentation\image\1.gif

Логистическая модель роста популяции

function LogModPop!(du,u,p,t)
    du[1] = p[1]*u[1]*(1-u[1]/p[2])
end
u0 = [39433556.0]
tspan = (0.0,100.0)
p = Float64[1.09, 15e7]
prob = ODEProblem(LogModPop!,u0,tspan,p)
sol = solve(prob,abstol=1e-6,reltol=1e-6, saveat=1.0)
R1 = [tu[1] for tu in sol.u]
plot(sol.t, R1, title="Модель Мальтуса", xaxis="Время",yaxis="Численность",label="u(t)", c=:blue)

anim = @animate for i in 1:length(sol.t)
    plot(sol.t[1:i], R1[1:i], title="Логистическая модель роста популяции", xaxis="Время",yaxis="Численность",label="u(t)", c=:blue)
end
gif(anim, "presentation//image//2.gif")

[ Info: Saved animation to D:\Education\КомпПрактикумПоСтатМоделированию\labs\gitrepo\lab6\presentation\image\2.gif

SIR

function SIR!(du,u,p,t)
    du[1] = -p[1]*u[1]*u[2] # S
    du[2] = p[1]*u[2]*u[1]-p[2]*u[2] # I
    du[3] = p[2]*u[2] # R
end
u0 = [39433553.0, 3.0, 0.0]
tspan = (0.0,20.0)
p = Float64[0.3,0.7]
prob = ODEProblem(SIR!,u0,tspan,p)
sol = solve(prob,abstol=1e-6,reltol=1e-6, saveat=0.01)
R1 = [tu[1] for tu in sol.u]
R2 = [tu[2] for tu in sol.u]
R3 = [tu[3] for tu in sol.u]
plot(sol.t, R1, title="SIR", xaxis="Время",yaxis="Численность",label="Susceptable", c=:green, leg=:topright)
plot!(sol.t, R2, title="SIR", label="Infected", c=:red, leg=:topright)
plot!(sol.t, R3, label="Recovered", c=:blue, leg=:topright)

anim = @animate for i in 1:length(sol.t)
    plot(sol.t[1:i], R1[1:i], title="SIR", xaxis="Время",yaxis="Численность",label="Susceptable", c=:green, leg=:topright)
    plot!(sol.t[1:i], R2[1:i], title="SIR", label="Infected", c=:red, leg=:topright)
    plot!(sol.t[1:i], R3[1:i], label="Recovered", c=:blue, leg=:topright)
end
gif(anim, "presentation//image//3.gif")

[ Info: Saved animation to D:\Education\КомпПрактикумПоСтатМоделированию\labs\gitrepo\lab6\presentation\image\3.gif

SEIR

function SEIR!(du,u,p,t)
    betta, delta, gamma, N = p
    s, e, i, r = u
    du[1] = -betta / N * s * i
    du[2] = betta / N * s * i - delta * e
    du[3] = delta * e - gamma * i
    du[4] = gamma * i
end
u0 = [0.98, 0.02, 0.0, 0.0]
tspan = (0.0,200.0)
p = Float64[0.8,0.4,0.3,1.0]
prob = ODEProblem(SEIR!,u0,tspan,p)
sol = solve(prob,abstol=1e-6,reltol=1e-6, saveat=0.1)
R1 = [tu[1] for tu in sol.u]
R2 = [tu[2] for tu in sol.u]
R3 = [tu[3] for tu in sol.u]
R4 = [tu[4] for tu in sol.u]
plot(sol.t, R1, title="SEIR", xaxis="Время",yaxis="Численность",label="Susceptable", c=:green, leg=:topright)
plot!(sol.t, R2, label="Exposed", c=:orange, leg=:topright)
plot!(sol.t, R3, label="Infected", c=:red, leg=:topright)
plot!(sol.t, R4, label="Recovered", c=:blue, leg=:topright)

anim = @animate for i in 1:length(sol.t)
    plot(sol.t[1:i], R1[1:i], title="SEIR", xaxis="Время",yaxis="Численность",label="Susceptable", c=:green, leg=:topright)
    plot!(sol.t[1:i], R2[1:i], label="Exposed", c=:orange, leg=:topright)
    plot!(sol.t[1:i], R3[1:i], label="Infected", c=:red, leg=:topright)
    plot!(sol.t[1:i], R4[1:i], label="Recovered", c=:blue, leg=:topright)
end
gif(anim, "presentation//image//4.gif")

[ Info: Saved animation to D:\Education\КомпПрактикумПоСтатМоделированию\labs\gitrepo\lab6\presentation\image\4.gif

Лотки-Вольтерры

using NLsolve
# Аналитическое решение
function find_equilibrium(a, c, d)
    function system!(du, u)
        du[1] = a*u[1]*(1-u[1]) - u[1]*u[2]
        du[2] = -c*u[2] + d*u[1]*u[2]
    end

    initial_guess = [0.5, 0.5]
    result = nlsolve(system!, initial_guess)

    equilibrium_point = result.zero
    return equilibrium_point
end
# Численное решение
function LotkiVolterry(a, c, d, x1_0, x2_0, dt, num_steps)
    x1 = x1_0
    x2 = x2_0
    results = [(x1, x2)]

    for _ in 1:num_steps
        x1_new = x1 + dt * (a * x1 * (1 - x1) - x1 * x2)
        x2_new = x2 + dt * (-c * x2 + d * x1 * x2)
        x1, x2 = x1_new, x2_new
        push!(results, (x1, x2))
    end  
    
    return results
end

a = 2.0
c = 1.0
d = 5.0
x1_0 = 0.15
x2_0 = 0.25
dt = 0.01
num_steps = 10000

results = LotkiVolterry(a, c, d, x1_0, x2_0, dt, num_steps)
R1 = [x[1] for x in results]
R2 = [x[2] for x in results]

equilibrium = find_equilibrium(2,1,5)

plot(R1, R2, title="Лотки-Вольтерры (фазовый портрет)", leg=:topright)
scatter!([equilibrium[1]], [equilibrium[2]], color="red", label="Точка равновесия")

anim = @animate for i in 1:length(R1)
    plot(R1[1:i], R2[1:i], title="Лотки-Вольтерры (фазовый портрет)", leg=:topright)
    scatter!([equilibrium[1]], [equilibrium[2]], color="red", label="Точка равновесия")
end
gif(anim, "presentation//image//5.gif")

[ Info: Saved animation to D:\Education\КомпПрактикумПоСтатМоделированию\labs\gitrepo\lab6\presentation\image\5.gif

Конкурентные отношения

function KonkOtn!(du,u,p,t)
    du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = -p[1] * u[2] + p[2] * u[1] * u[2]
end
u0 = [16.0,7.0]
tspan = (0.0,10.0)
p = Float64[0.02, 0.04]
prob = ODEProblem(KonkOtn!,u0,tspan,p)
sol = solve(prob,abstol=1e-6,reltol=1e-6, saveat=0.1)
R1 = [tu[1] for tu in sol.u]
R2 = [tu[2] for tu in sol.u]
plot(sol.t, R1, title="Конкурентные отношения", xaxis="Время",yaxis="Значение",label="x", c=:green, leg=:topright)
plot!(sol.t, R2, label="y", c=:orange, leg=:topright)

anim = @animate for i in 1:length(R1)
    plot(sol.t[1:i], R1[1:i], title="Конкурентные отношения", xaxis="Время",yaxis="Значение",label="x", c=:green, leg=:topright)
    plot!(sol.t[1:i], R2[1:i], label="y", c=:orange, leg=:topright)
end
gif(anim, "presentation//image//6.gif")

[ Info: Saved animation to D:\Education\КомпПрактикумПоСтатМоделированию\labs\gitrepo\lab6\presentation\image\6.gif

plot(R1, R2, title="Конкурентные отношения (фазовый портрет)", leg=:topright)

Консервативный гармонический осциллятор

function KGO!(du,u,p,t)
    du[1] = u[2]
    du[2] = -p[1]^2 * u[1]
end
u0 = [0.0, 1.0]
tspan = (0.0,10.0)
p = Float64[3.0]
prob = ODEProblem(KGO!,u0,tspan,p)
sol = solve(prob,abstol=1e-6,reltol=1e-6, saveat=0.1)
R1 = [tu[1] for tu in sol.u]
R2 = [tu[2] for tu in sol.u]
plot(sol.t, R1, title="Консервативный гармонический осциллятор", xaxis="Время",yaxis="Значение",label="x(t)", c=:green, leg=:topright)
plot!(sol.t, R2, label="v(t)", c=:orange, leg=:topright)

anim = @animate for i in 1:length(R1)
    plot(sol.t[1:i], R1[1:i], title="Консервативный гармонический осциллятор", xaxis="Время",yaxis="Значение",label="x(t)", c=:green, leg=:topright)
    plot!(sol.t[1:i], R2[1:i], label="v(t)", c=:orange, leg=:topright)
end
gif(anim, "presentation//image//7.gif")

[ Info: Saved animation to D:\Education\КомпПрактикумПоСтатМоделированию\labs\gitrepo\lab6\presentation\image\7.gif

plot(R1, R2, title="Консервативный гармонический осциллятор (фазовый портрет)", leg=:topright)

Гармонический осциллятор

function GO!(du,u,p,t)
    du[1] = u[2]
    du[2] = -2.0*p[2]*u[2] - p[1]^2*u[1]
end
u0 = [0.0, 1.0]
tspan = (0.0,10.0)
p = Float64[4.0, 0.2]
prob = ODEProblem(GO!,u0,tspan,p)
sol = solve(prob,abstol=1e-6,reltol=1e-6, saveat=0.1)
R1 = [tu[1] for tu in sol.u]
R2 = [tu[2] for tu in sol.u]
plot(sol.t, R1, title="Гармонический осциллятор", xaxis="Время",yaxis="Значение",label="x(t)", c=:green, leg=:topright)
plot!(sol.t, R2, label="v(t)", c=:orange, leg=:topright)

anim = @animate for i in 1:length(R1)
    plot(sol.t[1:i], R1[1:i], title="Гармонический осциллятор", xaxis="Время",yaxis="Значение",label="x(t)", c=:green, leg=:topright)
    plot!(sol.t[1:i], R2[1:i], label="v(t)", c=:orange, leg=:topright)
end
gif(anim, "presentation//image//8.gif")

[ Info: Saved animation to D:\Education\КомпПрактикумПоСтатМоделированию\labs\gitrepo\lab6\presentation\image\8.gif

plot(R1, R2, title="Гармонический осциллятор (фазовый портрет)", leg=:topright)